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Abstract 

We introduce a deterministic discrete-particle simulation approach, the Linearly-Transformed 
Particle-In-Cell (LTPIC) method, that employs linear deformations of the particles to reduce the 
noise traditionally associated with particle schemes. Formally, transforming the particles is jus- 
tified by local first order expansions of the characteristic flow in phase space. In practice the 
method amounts to using deformation matrices within the particle shape functions; these matri- 
ces are updated via local evaluations of the forward numerical flow. Because it is necessary to 
periodically remap the particles on a regular grid to avoid excessively deforming their shapes, 
the method can be seen as a development of Denavit's Forward Semi-Lagrangian (FSL) scheme 
[J. Denavit, /. Comp. Physics 9, 75 (1972)]. However, it has recently been established [M. Cam- 
pos Pinto, "Smooth particle methods without smoothing", arXiv:1112.1859 (2012)] that the un- 
derlying Linearly-Transformed Particle scheme converges for abstract transport problems, with 
no need to remap the particles; deforming the particles can thus be seen as a way to significantly 
lower the remapping frequency needed in the FSL schemes, and hence the associated numerical 
diff'usion. To couple the method with electrostatic field solvers, two specific charge deposition 
schemes are examined, and their performance compared with that of the standard deposition 
method. Finally, numerical Idlv simulations involving benchmark test cases and halo forma- 
tion in an initially mismatched thermal sheet beam demonstrate some advantages of our LTPIC 
scheme over the classical PIC and FSL methods. Benchmarked test cases also indicate that, 
for numerical choices involving similar computational eff'ort, the LTPIC method is capable of 
accuracy comparable to or exceeding that of state-of-the-art, high-resolution Vlasov schemes. 
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1. Introduction 



Although considered very efficient in many practical cases, particle-in-cell (PIC) simulations 
sometimes present levels of noise that make fine plasma phenomena very numerically expensive 
to resolve. The fact that particles are usually initialized with random procedures explains part of 
the statistical noise, yet there is another reason for the birth of strong oscillations in the numerical 
solutions. Indeed, it is known from the mathematical analysis of deterministic particle methods 
|l, 2] that a typical requirement for smooth convergence is that the radius e of the particles tend 
to at a much slower rate than the average grid spacing h used for their initialization, a property 
that is expensive to satisfy in practice. Here by "smooth convergence" we mean the pointwise 
convergence of the density function carried by the macro-particles, towards the exact solution / 
of the Vlasov equation. If the latter is a continuous function of the phase- space coordinates and 
if the convergence is pointwise, numerical solutions are indeed free of spurious oscillations, at 
least asymptotically. 

Specifically, smooth convergence requires s ^ with q < \, which can be interpreted as 
an extended overlapping condition: as the initialization grid gets finer, more and more particles 
must overlap. In PIC schemes (and more precisely, weighted PIC schemes with uniform Poisson- 
solver grids) the particle size is implicitly dictated by the J-dimensional mesh used in the field 
solver through its number of cells Nc ~ £~^, whereas the (average, if random) initial spacing 
can be derived from the number of particles A^p ~ k'^-d' ^ ^j^h d' denoting the dimension of the 
velocity variable. Therefore, to guarantee the smooth convergence of the numerical density one 
should increase the number of particles per cell consistent with the number of cells, i.e.. 



This exponent is always positive, and when d = d' ii is greater than unity, e.g., 1.5 if ^ = 
0.8. Hence for smooth particle simulations, the number of particles per cell should increase 
significantly faster than the number of cells. In practice such a condition is usually not met. 

On the mathematical level, particle methods that do not meet the extended overlapping con- 
dition may still converge towards a smooth / but only in a weak sense, i.e., in the sense that the 
local integrals of the particle density function tend to the same integrals of /. This case typically 
corresponds to simulations with strongly oscillating density functions, where accurate results 
can be obtained for certain integral quantities such as the electric field, or for the density itself, 
through appropriate smoothing procedures. And since the accuracy of the electric field is what 
matters most for the (electrostatic) dynamics of the system, strongly oscillating simulations can 
very well give satisfactory results on the longer scale sizes of physical interest. 

However, in cases where the physics of interest is in a region of low plasma density, smooth 
convergence seems to be necessary for precise measurements. For a variety of practical problems 
indeed (including backward Raman scattering plasma- wall transitions |0,|5|], halo formation 
in beams |6] and development of electron holes in the presence of a guide field f?]) physicists 
often need to resort to (grid based) Vlasov or Semi-Lagrangian solvers in order to obtain suffi- 
cient accuracy. Unfortunately these methods are known to be numerically expensive to run and 
challenging to implement, as they require the mesh to cover the whole phase space and can suff'er 
from diff'usive eff'ects. 



To reduce noise, Denavi t |[8|] p roposed a particle method later revisited as a Forward Semi- 
Lagrangian (FSL) scheme |9, 10], where the distribution function carried by the particles is 
periodically remapped to the nodes of a phase-space grid. This has a smoothing effect which 
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in practice eliminates the need for extended overlaping. However, frequent remappings can 
introduce unwanted numerical diffusion which in many cases contradicts the benefit of using low- 
diff'usion particle schemes. Resulting numerical diff'usion from the remappings can be reduced 
by use of high order adaptive schemes; see, e.g., |illlll2D. Other methods to reduce the noise 
have also been studied, such as wavelet-based denoising techniques; see e.g.. 

In this article, we present a new particle scheme where in addition to pushing the particle 
centers along their trajectories, one updates the particle shapes through the use of local linear 
transformations to better follow the local shear and rotation flows in phase space. As in the 
FSL scheme, the method is purely deterministic, and to prevent particles from being arbitrarily 
stretched, the particles need to be remapped periodically. However, significantly lower remap- 
ping frequencies are needed in practice, which results in higher accuracy and less numerical 
diff'usion. On a theoretical level this advantage is supported by the fact that for transport prob- 
lems with prescribed characteristic flow, the linearly transformed particle solutions are shown to 
converge in the uniform norm as h tends to 0, without any remappings; see 1 16]. 

Deforming the particles is not a new idea. For instance, our scheme can be viewed as a 
variation on Hou's formal vortex method 1 17] where the particles are deformed through a global 
mapping. In our method, each particle is transported by the linearized flow around its trajectory. 
Still on a formal level this approach coincides with a method presented by Cohen and Perthame 



111811 who established its first-order convergence, but they did not provide a numerical scheme 
to compute the deformation matrices. In the context of plasma simulation, an important class of 
deformed particle methods is off'ered by the Complex Particle Kinetic (CPK) schemes introduced 
by Bateson and Hewett In the CPK method, in addition to having the Gaussian shape of 

the particles transformed by the local shearing of the flow, the particles can also be fragmented 
to probe for emerging features, and merged where fine particles are no longer needed. When 
mature, our Linearly-Transformed PIC (LTPIC) scheme will incorporate some of the multilevel 
refinement features presented in [16], and the resulting adaptive scheme should be compared 
with the CPK method to determine which classes of problems best fit each method. 

The outline is as follows. In Section [2] we describe the LTPIC scheme for Idlv electrostatic 
plasmas: in Section [2Jl we introduce the main notations and present the general form of the 
numerical solutions. Although a wide range of shape-functions is supported by our approach, 
to illustrate the method we review in Section 12.21 one deterministic algorithm for initialization 
and remapping of B- spline particles, and we provide a correction scheme to make the remap- 
pings conservative. In Section [23l we then define a particle transport scheme that transforms the 
shapes in phase space, and is solely based on pointwise evaluations of the (forward) numerical 
flow. Two specific charge deposition schemes are then presented in Section [24l and a leap-frog 
time advance scheme implementing the method is described in Section 12.51 Numerical results 
involving several standard benchmark tests and a halo problem associated with an initially mis- 
matched thermal distribution are presented in Section [3l The results obtained demonstrate some 
advantages of our method compared to the classical PIC or FSL approaches. They also indicate 
that for numerical choices involving similar computational eff'ort, the LTPIC method is capable 
of accuracy comparable to or greater than state-of-the-art high-resolution Vlasov schemes. 

2. The numerical method 

To describe our method we may consider the normalized Idlv Vlasov-Poisson equation 



{dt + vdx -h E(x, t)dy} fit, z) = with r > 0, z = (x, v) g R^ 
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(1) 



which models the evolution of simplified plasmas and sheet beams; see for instance Ref. 1I21I1 . 
Here, x and v are dimensionless positions and velocities and £^ is a dimensionless electric field 
satisfying 

d^E{t,x)= I f{t,x,v)&v-ne, (2) 
where Ue is the density of a uniform neutralizing background cloud. 

2.1. Structure of the numerical solutions 

As in standard particle methods, we represent the phase- space density / with weighted col- 
lections of finite- size particles (index k) which are pushed along their trajectories correspond- 
ing to the discrete times tn = nAt, n = 0,1 ... ,Nt. However, in our method the particles also have 
their shape transformed to better represent the local shear and rotation flows in phase space, as 
illustrated in Figure [T] The particles can either be structured or unstructured. The first case cor- 
responds to the initialization and remapping steps, where particles are defined as tensor-product 
B- splines and centered on regular nodes 

^0 = (;^o^ ^0^^ = with G = {. . . , -1, 0, 1, . . .f. (3) 

Specifically, the univariate (i.e., one-dimensional) centered B- spline is recursively defined as 
the piecewise polynomial of degree p satisfying 

fl _i < r < i C'+l 
^o(x)^\' and Sp(x)^ ^ Sp.i(x)dx for p>l. (4) 

1^0 otherwise Jx--^ 

Thus Si(x) = max{l - |x|,0} is the traditional "hat- function", S3 is the well-known cubic B- 
spline supported on [-2,2], and so on, see e.g. Ref. 1I22I1 . The fundamental shape function is 
then defined on the two-dimensional phase space as a tensor product 

(p(z) = Bp{x)Bp{y) with support supp(^) = [-Cp,Cpf, Cp = (5) 

from which we derive a normalized (J(ph = 1), grid-scaled shape function (ph(z) = h~^(p(h~^z). 
Structured particles are then defined as translated versions of the latter, 

4,k(^) ^ <fh(z - zl) = h-^h-'z -k\ he 1}. (6) 

When transported by our method, particles become unstructured in the sense that their centers 
fl leave the nodes of the structured phase-space grid and their shapes are linearly transformed. 
That is, the positions of diff'erent parts of the "cloud" associated with a single particle advance 
with their own peculiar velocities, the velocities advance with their own peculiar accelerations, 
and the cloud distorts, but the distortion is constrained to be linear. Generic particles are then 
characterized by the 2 x 2 deformation matrices (initialized with Z)^ = ( q ? )) which determine 
the linear transformation of their shape, and numerical solutions take the form 

/;(^) = Z with <,(z) = ^,(Z)«(z-4)). (7) 

In Section [Z2I we shall describe a structured particle approximation operator 

H : f{z) ^ Yj ^k(f)n(z - zl) 
4 



acting on a generic density /, and in Section 12.31 and 12.51 we will construct a time-dependent 
transport operator 

The deformation matrices will be transformed with an area preserving scheme (det(Z)^^^) = 
det(Z)p = 1), so that the charges carried by the particles read (up to a constant factor) 

J wlifl,(z) dz = J wlifhiDliz - zD) dz = wl J ifhiz) dz = wl. 
In particular, the particle weights will not be modified by our transport operator. 




Figure 1: Structured particles (left) are defined at initialization and remapping steps, as tensor-product B-splines centered 
on regular nodes = hk, k e . Unstructured particles (right) are obtained by pushing the particle centers along their 
trajectories and transforming their shapes with a matrix representing the local Jacobian of the characteristic flow. 



To prevent the deformed particles from being arbitrarily stretched in one direction, we have 
chosen to periodically remap them onto the regular grid ©. Note that, as in Semi-Lagrangian 
methods, the remappings are likely to introduce unwanted numerical diff'usion. However, since 
our particle method is mathematically proven to converge without remappings |16], we expect 
the optimal remapping frequencies to be significantly lower than with Semi-Lagrangian schemes. 
This point will be numerically demonstrated in Section [3l 

The global structure of the scheme is then formulated as follows. First, a collection of 
weighted particles is initialized with 

^ = Ahfit = 0), 

then for n = 0, . . . , Afj - 1, we compute 

where /^J^^^^ if « > and mod (nA, A.,) = 0, 
■"' 1 /« otherwise 



to evolve the distribution with a given remapping period Ar^. 
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2.2. Initialization and remappings 

Since the fundamental shape function (ph is a B-spHne with the same scale as the spacing h 
of the regular grid ([3]), arbitrary polynomials with coordinate degree less or equal to p can be 
obtained by linear combinations of structured particles © derived by shifting iph on the grid, see, 
e.g., 1 22]. Therefore, to initialize and remap the particle densities we can use existing high order 
approximation schemes that rely on that property. One attractive method is given by the quasi- 
interpolation schemes described in Refs. |23] and |24]. Such schemes pass through data points 
when they are described by polynomial target functions / of a certain degree, and they have the 
advantage of computing high order B- spline approximants from local evaluations of the target 
function, unlike standard spline interpolation which requires solving a global system. Thus, in 
the univariate case the approximation aJ^^^^ takes the form 

A^l'^ : fix) ^ Yj ^k{f)^h{x - hk) with weights w,(/) ^ /z ^ ^/ fm + /)). (9) 

keZ \l\<mp 

Here we have denoted (ph{x) = h~^Sp(h~^x), and the ai = a-i are symmetric coefficients defined 
in such a way that A^^^^ f = f for any f(x) = ao -\- • • • -\- apX^. Specifically, they can be computed 
with the algorithm from Ref . ll23i Section 6] . For the first orders we find 

• nip = and ao = 1 f or = 1 , 

• mp = l and(ao,ai) = (|,-^) for = 3, 

• mp = 4 and (^0,^1,^2,^3,^4) = (ff, ^Is' sio' Tik) ^""^ P = ^^ 

In the bivariate case we can "tensorize" the above scheme, as it is easily checked that the operator 

Ah : f(z) ^ ^ Wk(f)(fh(z - Zk) with Wk(f) = ^ ^iK^l+i)^ = ^iM. (10) 

kel? \\l\U<mp 

reproducts any polynomial of coordinate degree less than or equal to p (here, ||/||co = max{|/;cL l^vl})- 
With standard arguments one can then show that the resulting approximation error converges as 
j^p+i smooth functions /, see e.g., Ref. [il6] . 

When the above scheme is used to remap a generic particle density of the form d?]), the 
mass of the resulting approximation is (setting a/ = for ||/||co > nip and using Yuiei? ^/ = 1) 

A,/«(z)dz = ^kiO r <Ph{z - zl)Az = ^kifD = h'Y = Y K' 

kel? ^ kel? k,l,k'EZ^ k'el? 

where w^, = h^w^, Z/gZ2 represents the charge deposited by the particle k' in the remap- 

ping process. Due to the shape transformations, this quantity generally diff'ers from the original 
w^, (indeed it vanishes if the support of misses the grid hi?). This shows that the quasi- 
interpolation is not conservative, but a locally conservative correction is easily implemented by 
depositing the local error w^, - w^, with a PIC-like method, which results in defining 

= ^kiO + Y « - K')h^vh{zl - 4,). (11) 

k'eZ^ 

Note that in practice the deposited fractions w^, can be evaluated by summing over / e 1} the 
values /./(^^) involved in the quasi-interpolation scheme. 
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2.3. Particle transport with linear transformations 

The LTPIC scheme is based on a LTP (Hnearly-transformed particle) transport operator Th[T] 
that transforms the particles through local linearizations of a given characteristic flow T. In this 
section we present the LTP transport operator in this general setting, and reserve for the following 
sections the description of the numerical flow. Schematically, one could indeed decompose the 
transport operator appearing in ([8]) as follows. 

1 . From the density carried by the particles one computes a numerical flow 

that approximates the (exact) characteristic flow of the Vlasov equation ([T]) over one time 
step [tn, tn+i\. Namely, the mapping : (x, v) i-^ (X, V){tn+i) that associates any phase- 
space point to the advanced-time point of its corresponding trajectory defined by 

{ X\t) = y(0 

{ with initial data (X, V){tn) = (x, v). (12) 

[ r(t) = E(t,X(t)) 

2. The particles are transported by the associated LTP transport operator. 

In Section [231 we will derive a leap-frog version of this approach, involving two intermediate 
flows r^"'^ and T"^'^ such that T"^'^ o^T^"'^ ^ T^^. The LTP operator wiU then be appHed twice 
per time step, i.e., we will define = ThVT^'^^ThVT^'^^. However, to simplify our presentation 
we shall consider in the remainder of this section that we are given a single flow ^ T^^. 

Applied to a generic particle ip\ ^ with deformation matrix the LTP transport operator is 



T^m - 'Pik^^hm--^)) ^ Kk=n(Dr(--zr)) with 



(13) 

where is a matrix representing the Jacobian of the flow at z^, defined as follows. An approxi- 
mated Jacobian matrix is first defined with a centered finite diff'erence scheme, 

(r,\j ^ (2h)-\(rj:)i{zl + hej) - {TJlUzl - hej)] ^ dj(r:M4) for 1 ^ ^ j ^ 2, (14) 

where we have denoted ej = (Sij)i<i<2. Here h is the grid spacing of the remapping grid, but 
a diff'erent spacing could be used as well. Next we observe that, while the exact flow has a 
Jacobian with uniform determinant equal to 1 , there is no reason why this should be true for the 
finite diff'erence approximation (fT4l) . To obtain a conservative transport scheme (in the sense that 
/ ^^;/(z) dz = f cpl/z) dz), we then define as 
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/^" = det(/p"/«. (15) 

To justify the above approximations inherent in Eqs. ([T3l)-([T5]) we temporarily assume that 
we can apply the exact flow 9^^^. Pushing a fixed- shape particle as in standard PIC schemes gives 

Tncir:^ : ^h(z - z^ ^ ^h(z - z^') with ^ r:^)- (16) 
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Now, since ([T2l) is reversible, the exact transport of an arbitrary phase-space density z) over 
the time step is f(tn+i,z) = f(tn, (^^x)~H^))- In particular, for the particle (fih(z - we have 

T,Ar:,] : niz - zD ^ (fh((r:j-\z) - zl\ 

to which ([T6b can be seen as the lowest order approximation. Enhanced accuracy is obtained by 
a first order expansion around z^: writing J<r{z) = {djTi{z))i<i y<2 the Jacobian matrix associated 
with an arbitrary flow : ^ M^, we let 

denote the linearized flow around z^. We then define a "formal" LTP transport operator as the 
exact tranport corresponding to this linearized flow, namely 

T^LTP [??x] = ^ex ] for the particle associated to (17) 

Applied to a structured particle, we observe that it reads 

W^Frx]:^/.fe-4)^^/.((/,Tkz-zr')) with zt'^r^M)^ Jk^Jmiziy (18) 

Now, replacing the flow 9^^^ by its numerical approximation 'T^" and using a finite diff'erence 
scheme for the forward Jacobian leads then to the practical LTP transport operator as defined 
by Eqs. (fT3l)-([T5l). A rigorous error analysis of this procedure is presented in |16]. This er- 
ror analysis demonstrates the global convergence of both the discrete ([T3l) and continuous ([TSl) 
schemes (without remappings) in the uniform norm, as h tends to 0, provided that the exact flow 
is approximated with sufficient accuracy. 

2.4. Conservative charge deposition schemes for linearly transformed particles 

To complete the LTPIC scheme we now describe how to compute the field from the linearly 
transported particles. For this purpose we equip the (Id) physical space with regular nodes 

Xi = iK with / G Z, 

and as in Section [Z2l we let ipw (•^) = ^^p(^) denote the scaled B-spline of mass 1 in the physical 
variable. Following the standard approach ll25ll we represent the charge density on the grid with 

p«(x)^^p«^,.(x-xO, (19) 

and solve the Poisson equation on the same grid. Specifically we represent the electric field with 

El{x)^Yj^'Eyh'{x-Xi) (20) 

with coefficients computed by a centered finite difference scheme, such as 

2h' ' {h'f 

Here the diff'erent normalization in ([T9b and (l2Qb allow the coeflftcients to be on the order of 
the point values E^^, (x/), whereas the coefficients are on the order of the local charges h' p\, (x/), 
consistent with standard notations, see I25j, Chapter 5]. To compute the local charges p^ several 
methods can be considered. 
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1. In the simplest approach the particles are seen as point particles and they deposit their 
charges in a way similar to PIC schemes, i.e., 

p-,^HYj<^h'{4-^i)' (21) 

Note that in this case the numerical scheme still differs from a standard PIC method, be- 
cause the particles are periodically remapped on the regular grid with a smoothing effect. 

2. To take into account the shape of the particles, specific deposition schemes can be used 
instead. They rely on an intermediate charge density defined as the exact integral of 
along the velocity variable, 

and on the use of univariate quasi-interpolation © to compute the local charges, 

Pl = 4'^~Pl i-e-, Pl^h'Yj aiPlM = ^' E Z ^I'PlkM- (22) 

\l\<mp kel? \l\<mp 

Note that a correction similar to (fTTI) can be used here to make the deposition conservative, 
in the sense that J p^(x) = J f^{z) dz. In the above deposition formula (l22l) , we observe 
that the evaluation of the "integrated particles" ^ is not straightforward, due to the linear 
transformation of their shape. To compute them we have considered two methods. 
2. a. In the first method (see Algorithm 12. II and Figure O we use a Gaussian quadrature 
p^^(x/+/) to evaluate each velocity integral p^^(x/+/). In order to be accurate this 
approximation requires a few quadrature intervals fitted to the particle support (pro- 
jected along the velocity variable) and a few Gauss points per interval. This makes it 
hard to apply in higher dimensions. 
2.b. In the second method (see Al gorithm 12.21 and FigureO we simply replace each inte- 
grated particle by a univariate weighted B- spline sharing the same first 3 moments. 
Specifically, we approximate p^ ^ as 

w'^ JC - jc" i 

pU^) - pU^) = ^) with 4, = h 7((D^)2,2)2 + ((D«),,2)2, (23) 

^h,k ^h,k 

so that ^ x^p^ j^(x) dx = x^p^ ^(x) dx holds for m = 0, 1, 2. The resulting imple- 
mentation is much simpler, and the extension to higher dimensions is straightforward. 

Algorithm 2.1 (Charge deposition with Gaussian quadrature). Let Ng and A^^ denote the pre- 
scribed number of quadrature intervals and Gauss points per interval in the v dimension, per 
particle. 

1. Loop over every active particle (p^ ^, i.e., over k ^I? such that 9^ 0. For conciseness we 
denote the deformation matrix by Z) = Z)^, and we observe that the particle support is 

supp(^^ ,) = (xl vD + D-\h[-Cp, cpf) = {(X, V) : \\D(x - 4, v - v^)|U < hcp}. (24) 
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2. Determine the x-projection [x^ - hx, -\- hx] of the support ([24b of ^, i.e., set 

hx = hxQiXn) = \ diam({4 + {D-^r)x : ||r|U < /ic^}) = /iC;,(|Z)2,2l + |/)i,2l), (25) 

where we have used Z)"^ = ( ^ since det(Z)) = 1. 

3. Loop over the non- vanishing point values ^(-^/O. namely over g Z such that 

14 -xv\< hx. i.e., r G { - /i,)J + 1, . . . , + /z,)] - 1 } . (26) 

Then define ^(x// ) ~ p^ ^(x/' ) with A/^g Gauss quadrature formulas using N'^ points in the 
V dimension: from (l24l) the interval [v~{V), v^{V)\ = supp(v i-^ (f^ ^(x//, v)) is given by 

vl + max,-.i,2{(A-,2)"HA-,i(4 - -^/O - ^^p) 



I v-(/0 = + max,.i,2{(A-,2)-HA-,i(4 - ^/O - /^c^p)} 
so that we may compute 



v^(n = + min,.i,2{(A,2)-HA,i(4 - ^/') + ^^/^)}' 



m'=0 m=l 



< f'\lk(^i'^v)dv=pl,(x,). 



(27) 



Here, Av = (v'^(r)-v (r))/NG, and Af, vf are the Gauss weights and nodes corresponding 
to the interval [0, 1], e.g., 

forA^^ = l: y^ = \, = 1, 

forA^^ = 2: yf = i(l ± = ^ /=1,2, 

r , .G _ 1 liG _ A 

^2 ~ 2' ^2 ~ 9 



- 2V 



Finally update the appropriate weights (initialized to 0) consistent with (l22l) , by setting 

p'l=p'l + Har-ipl^ixr) for i = T - I = T - m^, . . . , T ^ . (28) 

Algorithm 2.2 (Charge deposition with a moment method). 

1. Loop over the active particles (f^ ^, i.e., over k ^I? such that 9^ 0. 

2. Determine the support [x^ - /Zj^, -^^ + hx] of p^ ^, i.e., set hx = /Zx(^, ^) = ^^^/? with A^^ ^ 
defined as in (|23l) . 

3. Deposit the approximated charge contributions as above: loop over g Z satisfying 
and update the weights as in (l28l) , with the explicit expression (|23]) for p^ ^. 
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v+(/-2) 



v-(/-2) 




Figure 2: In the Gauss deposition scheme described in Algorithm 12. II the charge density p'lj^(x) = iR^hk^-^'^"^ 
associated to a linearly-transformed particle is deposited with a quasi-interpolation scheme where the required point 
values are estimated with Gaussian quadrature along v slices. 




Figure 3: In the moment deposition scheme described in Algorithm |2.2| no numerical integration is needed. Instead, the 
charge density pJJ ^(x) = ^ ^(x, v) dv is replaced by a B -spline pJJ^x) that shares its first 3 moments, and the latter 
is deposited with the quasi-interpolation scheme (compare with Figure [2). 



The good news is that in practice it does not seem necessary to resort to accurate piecewise 
Gauss quadratures. Indeed, in most of the numerical tests presented in Section [3] the results ob- 
tained with the simple moment deposition scheme (displayed) were compared with simulations 
using a Gauss deposition scheme with No = 4- quadrature intervals and N'^ = 3 Gauss points 
per intervals, and the differences were hardly visible. Maybe more surprisingly, we also com- 
pared these results with simulations using the much simpler, PIC-like point deposition scheme, 
and again the differences were hardly visible. This suggests that in many cases the oscillatory 
representation of the density (as seen by the field solver) still yields an electric field that is able to 
drive accurate dynamics. On one hand, we should not be overly surprised by such a fact, indeed 
it is routinely observed in PIC (and FSL) simulations. On the other hand, we could expect that 
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in some cases an accurate resolution of the field requires a smooth representation of the density 
as seen by the solver, especially when high order solvers are used. This is an important question 
that is familiar in the finite element community, and that shall be addressed in future research. 



2.5. Time marching scheme 

Equipped with the LTP transport operator Th[T] defined in ([T3l) for a generic flow and 
with the field solver described in Section 12.41 we are now in position to specify the numerical 
transport involved in the scheme ([8]), namely the operator 

J-h • Jh ^ h • 

To this end we consider a standard leap-frog time discretization. We first transport 

^mr;:'']/:' where r,''%,v)^(x+^v,v), 
then compute an intermediate electric field 

using finite diff'erences as described in Section [241 and finally complete the time step with 
/r ' = Tdr;:'']/:'' where r:'\x, V) = (x + f V, V = V + Af£«;kx)). 



3. Numerical simulations 

In this section we apply our LTPIC scheme on a series of Idlv test cases and compare the 
resulting solutions with classical PIC or FSL runs, the latter being obtained by freezing the 
particle shapes in our code - that is, by setting = ( q ^ ) for all k and n. To facilitate the 
comparison with LTPIC and FSL, we often indicate the number of particles used in a PIC run as 
a product (e.g., 128 x 128). This may correspond to a uniform grid where (weighted) particles 
are initialized, but in most cases the displayed PIC runs use unweighted particles initialized with 
a standard quiet start method. 

In Sections im and [T2I we first consider standard benchmark problems for which our results 



can be compared to the existing literature, see e.g. Refs. |26, 27, 28, 291 130|, |31|]. Next, in 
Section 13.31 we study a more applied test case consisting of a mismatched beam in a constant 
(continous) focusing channel, derived from the Id sheet beam model developed in Ref. |21]. 

For the benchmark test cases in Sections ITT] and [T2l we classically consider the normalized 
Vlasov equation ([T]) in x-periodic phase space [0, L] xR coupled with a periodic Poisson equation 



dMt,x)= I f(t,x,v)dv-ne, t > 0, xe[0,L]. (29) 

Here rie is the uniform (and constant) density of a neutralizing background cloud, and we com- 
plete ([29b with the standard condition ^ E(t, x) dx = 0. 
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3. 1 . Weak and Strong Landau damping 

We first consider the normaUzed Vlasov-Poisson system described above with perturbed ini- 
tial distribution ^ 

f{t = 0, X, v) = exp |-y j (l + A cos(y^jc)). (30) 

Consistent with classical benchmarks fl^ |27l Ell] we take k = 0.5 and set the perturbation 
amplitude A = 0.01 for the weak Landau damping test case, or A = 0.5 for the strong Landau 
damping test case (actually, for such a perturbation the field is only damped for times t ^ 10). 
In this section we use a cutoff' velocity Vmax = 6.5 and periodic boundary conditions at x = and 
X = L = In Ik. 

In Figure |4] we show the norms of the electric field (left panel, semi-log scale) and of the 
phase- space density (right panel) obtained with PIC and LTPIC simulations of the weak Landau 
damping. Results clearly show the noiseless aspect of the LTPIC method, as the theoretical 
damping rate (y = -0.1533) is matched with a low-resolution run using 64 Poisson cells and 
64 particles per cell. We also observe the classical recurrent relaxation occuring with period 
Tr ^ 60, in good agreement with the theoretical period L/Av ^ 62, see e.g., Ref. |26]. In 
contrast, a PIC run using the same number of cells and particles (labelled as PICi) is unable 
to predict the correct damping rate beyond ^ 5. And even with 1024 particles per cell (and 
significantly greater cpu time), the PIC3 run only predicts the correct rate until t = 20. We also 
see that the low-resolution LTPIC run does a significantly better job at preserving the norm 
of the density (a Vlasov invariant), compared to the low- and moderate-resolution PIC runs. For 
such a test case, only the high-resolution PIC run performs better with regard to the L2 measure. 

In Figure \5\ the same quantities are shown for the strong Landau "damping". Again, the 
low-resolution LTPIC run predicts the benchmarked rates for the initial damping and subsequent 
exponential growth. The low-resolution PICi run (with 64 cells and 64 x 64 particles) only 
predicts the initial damping. The moderate-resolution PIC2 run predicts correctly both rates 
(although with less accuracy for the growth period), but at a significantly higher cost in terms 
of memory and cpu time. As for the preservation of the norm of the density we observe that 
LTPIC does not perform better than PIC here, essentially due to the remappings. 

Finally, we found that the low-resolution FSL runs (using 64 cells and 64 x 64 particles) give 
energy curves very similar to the LTPIC ones, in both the weak and strong damping cases. These 
curves were omitted for readability. 

To better assess the noiseless aspect of our method, we also show in Figure[6]the phase-space 
density /^(•^, v) as it evolves in the time range tn e [0, 60], obtained with an LTPIC run with 
periodic (A^r = 4) remappings on a grid of 256 x 256 particles. Here the strong phase-space 
filamentation is accurately resolved. In particular it agrees very well with similar phase-space 
plots shown on Figure 10 in Ref. |[3lll . obtained with a high order Backward Semi-Lagrangian 
Discontinuous Galerkin (BSL-DG) scheme (our color scale is chosen in order to fit theirs). In 
that scheme the phase- space density is computed using a cartesian mesh of lower resolution in 
the X dimension (128 points), but same resolution in the v dimension (256 points), where the 
filaments are most difficult to resolve. Close examination reveals slightly better resolution of ffne 
structures with the BSL-DG scheme, but due to the local nature of the DG method, we note that 
for fffth-order accuracy each cell contains 15 basis functions, hence signiffcantly more (about 7.5 
as many) degrees of freedom are involved in the BSL-DG simulation. Comparison with another 
fifth-order BSL-DG simulation [30^ Fig. 7, bottom row] shows a better resolution for the LTPIC 
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20 40 60 80 100 120 20 40 60 80 100 120 

Time, Time, 

Figure 4: Weak Landau damping. L? norms of the electric field (left) and of the particle density (right) are plotted 
vs. tn = nAt for LTPIC and PIC simulations. All the runs use a time step of At = 1/8 and 64 cells for the Poisson solver. 
The PIC runs labelled as PICi, PIC2 and PIC3 use increasingly high numbers of particles, namely 64 x 64, 128 x 128 
and 256 x 256. The LTPIC run uses 64 x 64 particles and a remapping period A^r = 4. The approximate cpu times for 
these runs are 40 s (PICi), 90 s (PIC2), 330 s (PIC3) and 45 s (LTPIC). On the left panel the plotted slope (7 = -0.1533) 
matches the theoretical damping rate, see Refs. (^[3ll]. The quasi-periodic relaxation in the LTPIC curve is known as a 
Poincare recurrence, see text for details. 




10 20 30 40 50 60 10 20 30 40 50 60 

Time, Time, 

Figure 5: Strong Landau damping. Plotted quantities and numerical parameters for the PIC and LTPIC runs are the same 
than in Figure |4] The approximate cpu times for these runs are 20 s (PICi), 45 s (PIC2) and 24 s (LTPIC). The plotted 
slopes (71 = -0.2920 for the initial damping and 72 = 0.0815 for the growth between times t = 20 and t = 40) match 
benchmarked exponential rates, see e.g. Refs. 1271 , [3111 . 
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scheme, however in this case the BSL-DG run uses only half as many degrees of freedom as the 
LTPIC one. 

f[](x,v) at t, = f[^(x,v) at t, = 6 f[^(x,v) at t, = 9 




Figure 6: Strong Landau damping. Time evolution of the phase-space particle density fj^(x, v) obtained with an LTPIC 
simulation. This run uses a time step At = 1/8, a remapping period A^r = 4, a Poisson solver with 256 cells and 256x256 
particles. On the tn = 6 snapshot some oscillations are visible but they almost vanish in the subsequent plots, due to the 
particle remappings and the strong shearing of the flow. The approximate cpu times for this runs is 660 s. 

In Figures[7]and[8]we then compare how well different particle methods (namely PIC, LTPIC 
and FSL with various numerical parameters indicated in the figure caption) resolve the filaments 
in the strong Landau "damping" test case at ^ = 60. Phase- space densities obtained at ^ = 60 
with diff'erent schemes are shown on Figure[7l together with a reference solution obtained with an 
LTPIC run using improved numerical parameters relative to Figure [6l Again, the results clearly 
show that LTPIC and FSL are able to remove the noise. With PIC the localization of global 
patterns such as filaments and holes may be accurate, however the noise level is significant (and 
it remains so with finer simulations, not shown here). Results also show the eff'ect of varying 
the remapping period A^r in the FSL and LTPIC runs. For low remapping periods both methods 
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give similar results, which is expected since particles do not have time to deform much. For high 
remapping periods however the LTPIC performs significantly better: it introduces less diff'usion 
than FSL, and does not present the unphysical oscillations that start to appear in the filaments 
computed with the FSL method. This is also expected from the convergence analysis of the LTP 
transport operator ([T71) , which does not require remappings for asymptotic convergence. The 
good news is that this improved performance does not come at an expensive price: the measured 
cpu times are indeed similar for FSL and LTPIC runs, which indicates that the additional work 
of updating the deformation matrices does not represent a significant portion of the overall time. 



Reference solution unweighted PIC weighted PIC 




Figure 7: Strong Landau damping. Comparisons of phase-space densities obtained at tn = 60 with different methods. 
All the runs use a time step At = 1/8, a Poisson solver with 256 cells and 256 x 256 particles, except for the reference 
simulation, an LTPIC run with 512 cells and 512x512 particles. In the FSL and LTPIC runs the remapping period varies 
as indicated (in the reference run it is Atr = 4). The approximate cpu times for these runs are 4900 s (reference LTPIC), 
625 s (unweighted PIC), 650 s (weighted PIC), 690 to 720 s (FSL runs) and 655 to 665 s (LTPIC runs). 



Finally, in Figure[8]we show v-slices of the distribution at x = L/2 and tn = 60. Again, results 
show that the LTPIC scheme gives the best results: compared to the PIC method the noise has 
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been removed, and compared to the FSL scheme the numerical diffusion is significantly reduced. 



Reference solution unweighted PIC FSL with At,. = 4 LTPIC with At,. = 4 



Figure 8: Velocity profiles = L/2, v) at tn = 60 associated with some of the numerical solutions shown on Figure0 



3.2. Two-Stream instabilities 

Here we consider again the periodic Vlasov-Poisson system and set the initial distribution 
function as follows. 



1. (Weak instability.) First, to compare our results with Refs. |27l|280 we set 

2(1 + 5v') I / cos(2/:x) + cos(3/:x) 

f(t = 0, X, v) = -^^e 2 1 + A — + cos(^x) 

7V27r \ \ 1-2 

with k = ^ and a weak amplitude A = 0.01 for the perturbation. 

2. (Strong instability.) Next to compare our results with Refs. |l29ll3l|] we set 

fit = 0, X, v) = —=e 2 (1 - A cos(^jc)) 
V27r 

with k = ^ and a strong amplitude A = 0.5 for the perturbation. 

For the simulations we use a cutoff' velocity Vmax = 5 and periodic boundary conditions at x = 
and X = L = Injk. 

In Figure [9] we compare phase-space densities for the weak two-stream instability (case 1) 
obtained at r = 53 with PIC, FSL and LTPIC runs and various numerical parameters indicated in 
the figure caption. Again, the results lead to several observations. 

• First, LTPIC and FSL are able to remove the "noise" (i.e., the oscillations) for appropri- 
ate values of the remapping period A^r. In that regard our simulations show again the 
robustness of LTPIC compared to FSL, where strong oscillations appear for A^r ^ 2. 

• For low remapping periods FSL and LTPIC give similar results - an expected observation 
since particles have less time to deform. However a closer look at the filaments in the 
A^r = 1 case shows that the latter is less diff'usive. 

• Again, our measurements indicate that for similar numerical parameters the FSL and LT- 
PIC runs take similar computational time. This signifies that deforming the particles is not 
an expensive task in our code. 

• Finally we find that the LTPIC scheme is able to achieve the accuracy of some high- 
resolution state-of-the-art grid-based methods. For instance, the bottom left panel in Fig- 
ure [9] showing an LTPIC run using A^r = 1 and 128 x 128 particles is very similar to the 
right panel in Figure 11 from Ref. [28], obtained with a conservative third order WENO 
BSL scheme using a 256 x 512 phase-space mesh. 
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Reference solution unweighted PIC weighted PIC 




Figure 9: Weak two-stream instability (case 1). Comparisons of phase-space densities obtained at tn = 53 with different 
methods. All the runs use a time step At ^ 1/5, a Poisson solver with 64 cells and 128 x 128 particles, except for the 
reference simulation, an LTPIC run with 128 cells and 512 x 512 particles. In the FSL and LTPIC runs the remapping 
period varies as indicated (in the reference run it is A^r ^ 2). The approximate cpu times for these runs are 875 s 
(reference LTPIC), 30 s (unweighted PIC), 26 s (weighted PIC), 31 to 37 s (FSL runs) and 35 to 40 s (LTPIC runs). 



Turning next to the strong two-stream instability (case 2), we show in Figure [TOl the time 
evolution of the phase- space density obtained with an LTPIC run using 256 x 256 particles. Fine 
phase- space detail is resolved as the strong amplitude of the initial perturbation leads to filamen- 
tations. Again we can compare our results with high order state-of-the-art grid-based methods. 
For instance we observe that the bottom right panel in Figure [TOl showing the LTPIC density at 
t = 45 resolves the filaments with similar accuracy to that of the (center and bottom) panels in 
Figure 4 from Ref. |31], obtained by a fifth-order BSL-DG scheme using 129x 129 and 255 x 255 
phase-space cells, respectively. Here some complementary observations are in order. On the one 
hand indeed, a closer look at the bottom right panel shows that the LTPIC solution presents mod- 
erate oscillations in the inner filaments. Therefore, it is not strictly as accurate as the mentionned 
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BSL-DG simulations. (To remove these oscillations one may lower the remapping period A^r but 
at the cost of more diffusion; one should then use remapping operators less diffusive than the cu- 
bic spline quasi-interpolation.) On the other hand, as a forward particle method LTPIC is simpler 
to implement, and potentially cheaper to run compared to grid-based or BSL methods. Finally, 
we note once again that the LTPIC simulation shown in Fi sure [TOl involves significantly smaller 
numerical systems: as each fifth-order DG cell contains 15 basis functions in 2d, BSL-DG runs 
using 129 X 129 and 255 x 255 meshes involve approximatively 4 and 15 times as many degrees 
of freedom than the plotted LTPIC simulation. 



f[^(x,v) at t, = f[](x,v) at t, = 3 f[^(x,v) at t, = 6 




Figure 10: Strong two-stream instability (case 2). Time evolution of the phase-space particle density /^(x, v) obtained 
with an LTPIC simulation. This run uses a time step = 1/5, a remapping period A^r ^ 2, a Poisson solver with 64 
cells and 256 x 256 particles. The approximate cpu time for this run is 115 s. 

Therefore, in the top row of Figure [TT] we show PIC, FSL and LTPIC runs using 512x512 
particles. According to the above discussion, these high-resolution runs involve as many de gree s 
of freedom as the 129 x 129 BSL-DG solution shown in Figure 4, center row, from Ref. ll3in . 
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With such resolution, we find that the LTPIC achieves a higher level of details probably caused 
by reduced diff'usion, similar to the 256 x 256 BSL-DG solution shown in the bottom row of the 
same figure in Ref. iIsHi . In the bottom row of Figure [TT] we then show PIC, FSL and LTPIC 
runs using 256 x 256 particles, for comparison. To highlight the robustness of LTPIC compared 
to FSL the remapping period is taken higher than in Figure \W\ Results indeed show strong 
ripples in the FSL solutions, but almost none in the LTPIC ones. By running the FSL method 
with A^r = 1.8, we obtain solutions (not shown here) where oscillations are either significantly 
reduced in the 256 x 256 case or fully smoothed out in the 512 x 512 case. Finally, our cpu time 
measurements show that FSL and PIC runs using similar numbers of Poisson cells and particles 
require very similar computational eff'ort. With the same numerical parameters the LTPIC runs 
are only slightly longer, which again indicates that the extra work of deforming the particles is 
not dominant. 



51 2x51 2 PIC 512x512 FSL with Atr = 3 512 x 512 LTPIC with At^ = 3 




Figure 11: Strong two-stream instability (case 2). Comparisons of phase-space densities obtained at tn = 45 with high 
(top) and moderate (bottom) resolution runs. All the runs use a time step At = 1/5. The top runs use 128 cells and 
512 X 512 particles, with approximate cpu times 600 s (PIC), 630 s (FSL) and 680 s (LTPIC). The bottom runs use 
64 cells and 256 x 256 particles, with approximate cpu times 90 s (PIC), 90 s (FSL) and 110 s (LTPIC). Note that the 
remapping period At^ = 3 taken here for the FSL and LTPIC runs is slightly longer than in Figure [TOl 



3.3. Halo formation in a mismatched thermal sheet beam 

We now consider the case of a ID sheet beam in a continuous focusing channel with pre- 
scribed focusing strength 

<s)^kl, (31) 

as studied in Ref. |21]. Here no electron cloud is present (a: takes the role of a neutralizing 
species) and the density / = f{s, x, x') models an axially thin, transverse slice of a continuous 
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{djdz = d/dy = 0) ion beam composed of single species particles of charge q and rest mass m. 
The sHce propagates with velocity jStc = const and relativistic gamma factor = (1 - 
along the axial (z) direction. Here, c is the speed of light in vacuum. The beam phase space is 
described by the spatial coordinate x and the angle that the particle trajectories make relative 
to the longitudinal axis, and the independent timelike coordinate s represents the axial coordinate 
of a reference particle of the beam (or of the slice being followed), measured along the design 
orbit (nominally the machine axis). In this model the Vlasov equation reads 

d dH d dH d ] ^ 



with Hamiltonian 



ds dx' dx dx dx' 



2 2 my]pY' 
Here, the electrostatic potential is the solution to the Poisson equation 



-^{s,x) = -^ni{s,x) (34) 

solved with free space boundary conditions -^(±oo) = ±^^/ to obtain 

E{s,x) ^ ~is,x) = J( J' ni(s,x)dx-^Ni) (35) 

where ni(s, x) = f(s, x, x') dx' is the ion density in configuration space and Ni = ni(s, x) dx 
is the integrated ion density, or total number of ions - a constant, as particles are neither cre- 
ated nor destroyed. Following the procedure described therein, we shall first review how thermal 
equilibrium solutions can be obtained with physical scales roughly consistent with a recent exper- 
iment for beam driven Warm Dense Matter called the NDCX-I at Lawrence Berkeley National 
Laboratory |[32ll . Specifically, we shall consider a 100 KeV kinetic energy potassium K'^ ion 
beam (the axial velocity of which can be set nonrelativistically by mybl3lc^/2 ^ 10^ q) with the 
sheet beam perveance 

P-^^. (36) 

to be specified below. Following the analysis carried out in f2?], thermal equilibrium distribu- 
tions can then be obtained as follows. Given a specific value for the positive, dimensionless 
parameter 

,37, 

where cbp = [q^h/(€om)Y^^ is the plasma frequency formed from the peak density scale h, a 
normalized eff'ective potential is defined as the solution of the transformed Poisson equation 

i/^'^(x) = 1 -h A - e-"^^^^^ with boundary conditions i/^'^(0) = jAa(O) = 0. (38) 

Then, the thermal distribution given by 
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yields a stationary equilibrium solution f(s, x,x') = f^^(x, x') of the Vlasov-Poisson system ([32b , 
([33])-(l35l). Here, T* = T I {mybPlc^) is the dimensionless temperature associated with the thermo- 
dynamic temperature T (expressed in energy units), and Ad = [T /(moj^)]^^^ is the corresponding 
Debye length. We observe that the parameters in (|39b can be derived by first inverting (l37l) , i.e., 

n = 



^2(1 + A) ' 

and next infering from (l36l) and Ni = ni(x) dx = 2h ytAo ^"'Aa W t^^t 

, p_a^A)^ 



The resulting temperature is then 



py 1 + A 



The tune depression cr/cro - defined as the ratio between the phase advance of the particles 
in the presence and absence of beam charge - can be calculated |21] as 



cr/cro 



1 



1 (f^^e-^^^^Hx) 



3/2 



V3(l+A)(J^"jc2^-^A(x) 



(41) 



By solving numerically (1381) , (|4T1) , it is then possible to prescribe a specific tune depression 
and derive the corresponding value of A to specify the equilibrium distribution: the resulting 
parameters are given in |21, Table II] for regularly spaced values of cr/o-Q e {0.1, 0.2, . . . , 0.9}, 
and in Table[T]below for a strong tune depression, i.e. cr/o-Q = 0.1. For the purpose of comparing 
our results to typical NDCX-I experiments, we set the focusing strength to k^^ = cro/Lp in such 
a way that free particles have a phase advance of ctq ~ n/3 per lattice period Lp = 0.5 m, and 
we set the perveance by taking P/k^^ = 0.01. We note that for a sheet beam, the perveance has 
dimension 1/length and P/k^^ is dimensionless. The distribution corresponding to cr/o-Q = 0.1 
corresponds to a highly nonlinear form in x due to the radial beam extent and the nonlinear 
solution for the eff'ective potential j/^a- 



depression 


parameter 


temperature 


Debye length 


rms radius 


peak density 


cr/cro 


A 




Ad 




h 


0.1 


5.522 X 10"' 


3.463 X 10"^ 


2.810 X lO-^m 


4.822 X 10-^m 


4.848 X 10^^ 



Table 1: Physical parameters for the matched thermal sheet beams with 100 KeV ions, corresponding to an axial 
beam velocity of /3bC with /3b = 2.343 x 10~^ and relativistic factor ~ 1. We take k^Q ^ 2.094 corresponding to a 60° 
phase advance per lattice period Lp = 0.5 m, and the perveance is set to P = O.OI^^q . 

Finally, from (|39l ) we derive an initially "mismatched" beam through a canonical transfor- 
mation that dilates the distribution in the spatial dimension while preserving its perveance and 
initial eff'ective phase- space area (emittance) by taking 

f(s = 0, X, /) ^ 1-,/ux] , > 0. (42) 



Here // corresponds to the mismatch parameter, defined as the ratio of the initial (rms) beam 
radius to the radius of the matched beam, see e.g. Ref. |33]. 

In Figure [12] we show the evolution of a mismatched beam with a thermal equilibrium form 
specified by (l39l ) using the procedure outlined above with jd = 1.25 and a tune depression of 
cr/o-Q = 0.1. Here the numerical solution is computed with an LTPIC scheme on a computational 
domain corresponding to \x\ < 15mm and \x'\ < 14.5mrad. In the phase-space plots we vizualise 
the tenuous halo that evolves from the initial distribution by taking contours of the numerical 
density using exponential increments. Filled color contours illustrate the core of the phase- space 
density. 

In Figure [13] we next compare the phase- space densities at s = 20m with halo contours 
obtained with diff'erent schemes using 256 x 256 particles. Here we observe that the unweighted 
PIC has a low level of noise in the core but misses almost all of the halo. The weighted PIC 
simulation catches a fair proportion of the halo but in a very fragmented way, and in addition it 
has a high level of noise in the core. In contrast, the FSL simulation with short remapping periods 
(shown on the center left panel) does a reasonable job, although it still misses some part of the 
halo arms. For longer remapping periods (center and center right panels) it is severly hampered 
by phase-space oscillations. In the LTPIC simulations (bottom panels) these numerical artifacts 
are significantly reduced, which shows once more the ability of this new approach to remove the 
noise at reasonable computational cost, and with similar or improved accuracy. 

4. Conclusion 

We have presented a new deterministic PIC method for electrostatic plasma simulations, 
wherein finite- sized particles have a shape function that is linearly transformed in time to ap- 
proximately follow the flow and thereby reduce the oscillations traditionally observed in standard 
PIC simulations. Although this method may be seen as an extension of the remapped-particle 
(FSL) scheme |8] (due to the practical need to remap the particles after the flow has evolved 
significantly), its relative robustness to low remapping frequencies makes it actually closer to a 
proper particle scheme. By testing our method on benchmarked test cases we have demonstrated 
its ability to eff'ectively reduce the noise and reach accuracy levels similar to those of expensive 
high order state-of-the-art Vlasov schemes. 

Simulations of real- world systems often need to capture the production of phase space struc- 
tures by collective interactions. For example, in accelerator physics, beam halo must be min- 
imized in order to limit the unwanted particle loss on machine surfaces. The simulated halo 
location and density must be quantitatively correct, in a calculation wherein the far denser core 
of the distribution must also be evolved self-consistently (since the core fields influence the halo 
dynamics). On test problems of this type the LTPIC method performs well, and we believe that it 
has considerable promise to augment the standard PIC methods predominately employed today. 
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Figure 12: Mismatched thermal sheet beam in a continuous focusing lattice with phase advance of 60° per period 
Lp = 0.5m, and tune depression (t/o-q = 0.1. The plots show the evolution of the phase-space density (with respect 
of the timelike, longitudinal coordinate s) obtained with an LTPIC simulation. This run uses a time step = Lp/ 16, 
a remapping period A^r = 2.5Lp, a Poisson solver with 128 cells and 256 x 256 particles. The halo is shown through 
isolines corresponding to values of 10~^, 10~^, . . . , 10~^ of the peak initial density. The approximate cpu time for this 
run is 220 s. 
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Figure 13: Mismatched thermal beam. Comparisons of phase-space densities obtained at Sn = 40Lp with different 
methods. All the runs use a time step As = Lp/ 16, a Poisson solver with 128 cells and 256 x 256 particles, except for the 
reference simulation, an LTPIC run with 256 cells and 512 x 512 particles. In the FSL and LTPIC runs the remapping 
period varies as indicated (in the reference run it is A^r = 2.5Lp). The approximate cpu times for these runs are 1450 s 
(reference LTPIC), 225 s (unweighted PIC), 145 s (weighted PIC), 220 to 265 s (FSL runs) and 200 to 245 s (LTPIC 
runs). 
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